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SUMMARY 


This report presents a new approach to the solution of matrix problems 
resulting from integral equations of mathematical physics- Based on the inherent 
smoothness in such equations, the problem is reformulated using a set of ortho- 
gonal basis vectors leading to an equivalent coefficient problem which can be of 
lower order without significantly impairing the accuracy of the solution. This 
approach has been evaluated using a two-dimensional Neumann problem describing 
the inviscid, incompressible flow over an airfoil. Two different kinds of mode 
functions have been investigated, namely eigenvector series and Fourier series. 

The first method has the disadvantage that the calculation of each eigenvector 
is a time-consuming process, and the results obtained do not indicate that this 
is offset by the need to use only a small proportion of the eigenvectors. On the 
other hand, the Fourier method is attractive since all of the expansion coefficients 
can be calculated very rapidly using the Fast Fourier Transform algorithm. The 
method developed here uses all of these coefficients in an approximate method 
which exploits the known structure of the transformed coefficient matrix and very 
promising results for the flow over a realistic airfoil are obtained. The mode- 
function method is aimed at reducing the computer time for the solution of the 
equations for large three-dimensional cases. On the basis of the results pre- 
sented here, it is shown that an order of magnitude reduction in this computer 
time can be expected for such problems as compared with the time for a direct 
matrix solution. 


V 


T.O INTRODUCTION 


1.1 General Background 

The result of the discretization of the governing integral equation for 
many physical problems is a large system of linear simultaneous equations. Such 
systems are usually highly "nonrandom" since the original physical problem gen- 
erally possesses a smoothly varying bahavior away from isolated singular regions. 
This is particularly true of the integral equations of many aerodynamic flow 
problems where the complexity comes from the geometry and a large number of 
points must be used to represent this adequately. This large number of points 
leads to a smoothly varying local behavior of the discrete quantities, but it 
also leads to the high computing costs associated with such calculations. For 

m simultaneous equations, the computing time required for a direct solution is 

3 2 

proportional to m /3 while that required for an iterative method is pm 

where p is the number of iterations required to obtain a fully converged solu- 
tion. On the other hand, the central thesis of the mode-function approach pre- 
sented here is that, for problems of this kind, advantage can be taken of the 
locally smooth character by replacing such a matrix equation with a related coef- 
ficient problem for certain preselected mode functions. Since a smooth function 
can, in general, be represented adequately in terms of a relatively small number 
of such functions, the size of the coefficient matrix can be significantly 
smaller than that of the original problem without seriously affecting the 
accuracy. 

One problem which is of great practical significance is the calculation of 
the internal or external, inviscid, incompressible flow about a given shape. This 
is the classical boundary- value problem for Laplace's equation with a well- 
established mathematical pedigree (the Neumann problem) which leads to an integral 
equation of the second kind. "Panel methods" which provide a numerical solution 
of this problem have been used for a number of years for aircraft design both in 
the United States and in Europe. The main difficulty associated with their use 
is the cost incurred both in terms of the labor involved in the preparation of 
the input data and the computational effort involved in the solution of the 
resulting equations. The reduction of this computation time is, therefore, a 



primary aim of the mode-function concept although it should be emphasized that 
the method presented here is applicable to the solution of any matrix problem, 
provided that the matrix possesses a dominant diagonal. The panel method 
developed at Douglas Aircraft Company is based on a source distribution and the 
associated influence matrix does possess a dominant diagonal. This problem is, 
therefore, particularly suited to the mode -function approach and so it has been 
used as the basis for most of the results presented here. 

A study of the application of the mode-function method to two-dimensional 
Neumann problems has been undertaken although the main practical application of 
the techniques would be in three-dimensional calculations for which the computa- 
tion times are significant. However, this present study demonstrates the valid- 
ity of the mode-function concept while providing a guide to the reductions in 
computing time which can be expected for three-dimensional calculations. This 
work is also highly relevant to three-dimensional problems since, as suggested 
in section 6, the fitting of a large matrix associated with a complex three- 
dimensional configuration would be undertaken in smaller blocks corresponding 
to two-dimensional sections. The operation count presented in section 6.2 
indicates that a three-dimensional mode-function method could offer significant 
savings in computer time compared with a direct matrix solution method. 

Before proceeding with an outline of the theoretical basis of the mode-func- 
tion approach, it will, however, be worthwhile to discuss the principal features 
of the panel method developed at Douglas Aircraft Company. 

1.2 The Hess Panel Method 

The low-speed flight regime occupies a special place in the field of aero- 
dynamics; first because it is the one regime in which all vehicles must operate, 
and second because only at low speed is the governing flow equation linear 
(Laplace's equation) so that powerful flow-calculation techniques can be employed. 
In particular, the problem can be formulated as an integral equation for a certain 
singularity distribution over the surface of the body about which flow is to be 
computed (ref. 1). This procedure owes its great efficiency to the fact that the 
domain of calculation can be restricted to the body surface, i.e., calculations 
need not be performed in the field of flow. In three-dimensional problems the 
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numerical implementations of this procedure have represented the body surface by 
a large number of small four-sided surface elements or "panels" and thus have 
come to be called "panel methods." In the Hess version, these panels are placed 
on the actual body surface. This renders the method numerically "exact" and 
applicable to any arbitrary body. Some alternative formulations place panels 
interior to the body, e.g., on the mean surface of a wing. These latter methods 
are applicable only to a restricted class of bodies. 

The above-described method of approximating a three-dimensional body is 
shown in figure 1. On each element, a control point is selected where the bound- 
ary condition of zero normal velocity is to be applied and where surface veloc- 
ities are ultimately calculated. A "matrix of influence coefficients" is then 
calculated. This consists of the complete set of velocities induced by the ele- 
ments on each others* control points for unit values of all singularity strengths. 
The integral equation that expresses the zero normal- velocity boundary condition 
is then approximated by a set of linear algebraic equations for the values of the 
singularity strengths on the surface elements. In lifting cases, the Kutta con- 
dition along the trailing edge yields a small number of additional equations. 

The order of the "matrix of influence coefficients" and the order of the coef- 
ficient matrix of the linear equations are both equal to the number of elements 
used to approximate the body (or nearly equal, depending on details of the numerical 
procedure). In the basic form of the method, (ref. 2) the surface elements are 
plane quadrilateral "panels," and each singularity strength is assumed to be 
constant over each element — a step function from element to element. The two 
time-consuming parts of the computing task are the calculation of the "matrix 
of influence coefficients" and the solution of the resulting linear equations. 

The relative importance of these two tasks depends on the body considered, the 
element number and the type of numerical technique employed. However, as a 
general rule, for complicated bodies which require very large numbers of panels, 
the solution of the linear equations is the harder task. 

The Hess three-dimensional panel method in both its lifting (ref. 2) and non- 
lifting (ref. 3) versions has been distributed very widely to various government 
agencies, industrial concerns and universities. In general, very satisfactory 
results have been obtained. Probably the principal disadvantage of this or any 
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three-dimensional panel method is the expense of running each case. For a 
complete aircraft configuration, a very large number of panels is required to 
adequately represent the geometry, and it is exoensive both to generate the 
necessary geometric input data and to solve the resulting system of equations. 

A major reduction of the time for the preparation of the input data has been 
achieved by the automatic geometry handling programs developed by Halsey (ref. 4j. 
This report is concerned with the second problem, that of reducing the computina 
cost of the matrix solution. 

1.3 The Mode Function Approach 

Having outlined the panel method to which the new approach is directed, a 
more detailed presentation of the mode-function theory can now be given. The 
fundamental assumption on which it is based is that both the singularity strength 
and the elements of the influence matrix must necessarily be slowly varying func- 
tions of position on the body surface. Therefore, their values on adjacent 
panels cannot be very different, so that solving the complete set of linear equa- 
tions in which every unknown is assumed to be independent is wasteful in the 
sense that not all of the available information is being used. In the original 
problem, neighboring quantities are not in any sense independent. The mode-func- 
tion method, therefore, seeks to fit the matrix and the singularity distribution 
to account for their continuity and dependency. The new problem becomes that of 
determining the "strength coefficients" of the chosen mode functions. Since, 
for example, a relatively small number of Fourier coefficients will contain as 
much information as a larger number of data points, this coefficient problem can 
be of much lower order than the original problem. 

Two basic kinds of fit have been investigated based on Fourier series and 
on eigenfunction expansions. While the eigenfunction method appears to offer 
the advantage of using mode functions specifically related to the individual 
problem, they have the disadvantage that their calculation is a time-consuming 
process. Their use would, therefore, only be justified if accurate solutions 
could be obtained in terms of a very small number of functions and the results 
presented here indicate that this is not achieved in practice. On the other 
hand, the use of Fourier series is attractive since these functions are well 
understood and the use of the Fast Fourier Transform (FFT) algorithm offers an 
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extremely efficient way of calculating all of the coefficients for a given 
function. The details of the application of each of these methods will be 
presented in the next section and some results bearing out these statements 
will be presented in section 3, However, before proceeding to that 
there are a few general points on the application of the mode-function method to 
two- and three-dimensional Neumann problems which should be discussed. 

The basic integral equation for the source distribution on a three- 
dimensional nonlifting body takes the form (ref. 2) 

cr(p) - j^a(q)K(p,q)dS = -n(p) • (1) 

s 

where a(p) is the unknown source density at a point p on the body surface S. 
The kernel function K(p,q) represents the negative of the normal velocity 
induced at the point p by a unit source distribution over the surface element 
dS at the point q. 

The corresponding equation for a two-dimensional problem takes the closely 
related form 


cf(p) - y*cf(q)K(p,q)ds = -n(p) 


V 


( 2 ) 


where the kernel function K now takes on a simpler form. However, when equa- 
tions (1) and (2) are discretized, both lead to a matrix equation of the form 


m 


E 

j=i 




f. 

1 


i = 1 , ... m 


(3) 


where m is the total number of surface panels or arc length elements. The 
values A-j are known as the influence coefficient matrix. Thus, it can be 
seen that, in principle, the two- and three-dimensional problems are mathe- 
matically equivalent, although in general the value of m will be much larger 
and the calculation of the influence coefficients will be more complicated for 

three-dimensional flows. Also, the structure of the matrix [A..] and the 

^ 0 
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column vectors [a.] and [f,] will be rather different in the two cases. 

J * 

This point will, however, be considered in greater detail in section 6 when 
the extension of the mode-function concept to three-dimensional problems is 
discussed. The theory and results presented here will therefore be directed at 
the two-dimensional problem. 


The standard solution method based on a source distribution combines the 
"self-influence term" with the remaining terms to recast equation (3) as 





i = 1 , . . . m 


(4) 


where This equation now takes on a form more closely related 

to an integral equation of the first kind. However, since A.. << 1, it should 

^ J 

be noted that it still possesses the essential qualities associated with integ- 
ral equations of the second kind, namely that the matrix [A--] is diagonally 

• ij 

dominant. Therefore, it has desirable characteristics for solution both by 
direct Gaussian elimination and by iterative methods. It will be shown in 
section 4,5 that this property is also necessary for the success of the mode- 
function method presented here. 


From the mode -function 
values of a function of two 
a large spike corresponding 
matrix in terms of smoothly 
treat the formulation given 


point of view, the matrix A!., regarded as discrete 
variables, has the undesirable property that it has 
to the dominant diagonal. Since we seek to fit the 
varying mode functions, it is therefore natural to 
in (3) and attempt to fit the matrix A... 
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2.0 GENERAL THEORY 


2.1 The Eigenvector Form 

Having established the background for the mode- function concept, it is now 
appropriate to present the theory and establish a notation for a more detailed 
discussion. The eigenvector method will be presented first since this can pro- 
vide some insight into the question of whether to use a formulation based on 
equations of the first or the second kind. Despite the close relationship 
between such equations, it should be pointed out that this is not simply a ques- 
tion of whether or not the diagonal term is included in the matrix. If a normal- 
velocity boundary condition is used, a pure source distribution leads to an 
equation which has the natural form of an equation of the second kind, whereas a 
pure vortex distribution will lead to an equation of the first kind. 

Matrix notation will be used throughout since some of the more involved 
aspects of the Fourier version to be presented in subsequent sections can most 
easily be condensed into such a format. Thus, equations (3) or (4) can be 
represented as 

f = e_ — A£_=A’p^ (5) 

where f^ and are column vectors consisting of the given "onset flow" and the 
unknown singularity distribution while A and A' are the two different forms 
for the influence coefficient matrix. The precise form of these matrices will 
depend both on the type of singularity distribution selected and on the way in 
which the problem is discretized. However, to simplify the discussion, it 
will be expedient to introduce the base case in which p is the source density, 
assumed to be constant over each element, and A is obtained by integrating 
the kernel over linear elements. The modifications introduced by the applica- 
tion of the mode-function method to alternative singularity distributions, pre- 
sented in section 5, is straightforward. 

Pursuing first the idea of fitting the matrix A in equation (5), note that 
Tt can be regarded as a function of two variables defined at m equally spaced 
intervals, where the independent variables correspond to the row and column 
numbers. By separation of variables, this function can be approximated by a 
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product of two functions £= (f-j , f2» ... f^) and 3^* (g^ , 92* ... g^^) 

dependent on the row and column number, respectively. Thus, an element A., 
of the A matrix is approximated by f.g. and the square of the error involved 

* U 

in this approximation is 

m m 
j=l i=l 

Differentiating this expression with respect to f^- and g^ leads to the con- 
ditions under which this error will be a minimum. Thus 


Aa = k^f , 
and 

Af = h\, 


2 2 2 2 

= f? + f? + . . . + 

I c m 


( 7 ) 


where A is the transpose of A. Therefore, f^ and £ will be eigenvectors 
of the matrices AA and AA since 


AAf = h^k^f 

and (8) 

AAa= h^k\ 

Now AA and AA are syimietric matrices and so, provided that their eigenvalues 
are distinct, they will each possess a set of orthogonal eigenvectors which must 


satisfy 

''2k = ^kik 


and 


( 9 ) 


*fk ' 'k^k 


or 

- ik 


and 




, k = 1, 2, ... m (10) 


The corresponding fit in terms of the first n of these eigenvectors is 

A = + X2f2ig + ... + 
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There are several points about this scheme which should be made at this stage. 
The procedure would be out of the question if the eigenvectors had to be cal- 
culated directly from equation (10) since the formation of the matrix products 

~ o 

AA and AA would require 0(m ) operations, and so this operation alone would 
require more computer time than the direct solution. However, it is possible to 
show that an adaptation of the traditional power method of iteration can be 
applied to equation (9), and one eigenvector can be calculated in 0(pm ) oper- 
ations where p is the number of iterations required (typically p '^^5 or 6). 
It is therefore clear that, in order for this method to be competitive with the 
direct matrix solution, n must be small compared with m since the operation 
count could otherwise approach or exceed 0(m /3) which is required for the 
direct solution. 


In terms of the matrix notation, equation (11) can be written 

A = FDG (12) 

where D is the diagonal (n x n) matrix of eigenvalues and F and G are 
the basic (m x n) “fit matrices" whose columns are defined by the eigenvectors 
f|^ and 3 ^ (k = l,...n), respectively. Similarly, one can express the 
vectors £ and in terms of the matrix G to give 


f_ = G^ and 

so that equation (5) can be written 

G^ = G^ — FDGGa^ 


(13) 

(14) 


It has been noted that the eigenvectors f^ and form two orthogonal sets 
of vectors and so, choosing them to have unit length, it follows that the matrix 
G is orthogonal, so that GG = I. Equation (14), therefore, gives 


^ where C = GFD (15) 

3 

The solution of this equation will involve 0(n /3) operations, but we have 
already asserted that the scheme will only be competitive if an adequate approxi- 
mation to equation (5) can be obtained for n << m. Thus, the solution time 
for equation (l5) will be negligible compared with the time required to find 
the eigenvectors. 


The accuracy of this truncated equation will in turn be governed by the fit- 
tability of the matrix A and the nature of the vectors f and £_. Since the 
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behavior of these eigenvector expansions is unknown a priori, the validity of thi 
truncation must be based primarily on the numerical results obtained. This situ- 
ation contrasts strongly with Fourier series where the behavior and the factors 
governing convergence are well understood. 


Before proceeding with an outline of the theory behind the Fourier series 
expansion, it will be worthwhile to derive the truncated eigenvector equation 
corresponding to an integral equation of the first kind since this throws some 
light on the general fittability of such equations. If F' and G‘ are the 
eigenvector matrices corresponding to the matrix A' then, following equations 
(12) and (13) 

A' = F'D'G' 
and 


f = F'd‘ and p_ = 

Equation (5) can then be written as 


or 


d‘ = D'a' 




k = 1 , ... n 


(16) 


(17) 


It follows that, since the eigenvalues decrease as k increases, the 
expansion for p will only converge if the numbers dj^ decrease even more 
rapidly. The form of equation (17) does indicate, however, that for problems 
which fall naturally into the form of an equation of the first kind, the converg- 
ence of the eigenvector expansion of the solution will be poor. This is con- 
sistent with the results derived from the Fourier expansion method for which 
unsatisfactory results were obtained when applied to an equation in this form. 

In fact, the method presented in section 4.5 is shown to be valid only for equa- 
tions of the second kind for which the corresponding matrix A' possesses a 
dominant diagonal. 


2.2 The Fourier Form 

Much of the discussion given in the previous section is also applicable to 
the Fourier format, and so the theory can be presented more briefly. It should 
be noted at the outset that this mode of fitting relies on the existence of 
known rapid calculation methods and known artifices to accelerate converaence 
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even though Fourier series as such have no particular relationship to Neumann 
kernel functions. It is, however, based on the observation that the right-hand 
side vector f_, the source vector g_ and the influence matrix are arrays of 
numerical values which are necessarily periodic since f^ = ^m+1 * ^1 " '^m+l ’ 


N,1 ' ^i.m+1 


and A-j ^ = A^^-j y For instance, in figure 2, with m = 12 the 


points i = 1 and i = 13 will be identical. The analysis presented in this 
section will be based on this periodic behavior although it should be noted that 
in practice there may arise discontinuities in the fitted vectors such as that 
occurring across the sharp trailing edge of an airfoil. Modifications to the 
fitting procedure necessary to cope with such discontinuities will be presented 
in section 4. 


The best rate of convergence for a finite Fourier series is obtained by 
using a pure sine series. However, such a series will vanish at the ends of 
the interval over which the series is defined, and so to fit a general function 
an additional constant must be introduced. Thus, if the function to be fitted 
is defined by the values (f-j, f2» ... f^) the series expansion will be 

m 

fj. = b^ + ^ ^ bt^ sin(j'l ) (k-1 )h where h = ir/m (18) 

k==2 

In this expression, b^ = f^ and the remaining coefficients are the standard 
Fourier sine coefficients of the function values (f^ — b^). Using the matrix 
notation introduced in the previous section, this can be written 

f = Gb (19) 

in which the transformation matrix now has the form 


G = 


1 

1 

1 


’ll 


’21 


’12 

’22 


0 

^1 ,m-l 
^2,m-l 


(20) 


^ Vl,l Vl,2 


S 


m-1 ,m-l 
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where S-. = sin(ijh), h = 7r/m. Similarly, the source density can be 
expressed as a sine expansion with coefficients 

p = Ga (21 ) 

Application of this process first to the rows of A and then to the columns of 
the resulting matrix leads to 

A = GBG (22) 

where B is a matrix consisting primarily of Fourier coefficients. 

The elements of ^ and B can be obtained by the application of the 
FFT algorithm to one- or two-dimensional arrays. Substitution of these expres- 
sions into equation (5) gives 

= G^ — GBGGa (23) 

Let D - GG, then D can be shown to have the form 



where 

- cot(kh/2) (k odd) 

=0 (k even) 

The simple structure of the matrix D will make the calculation of the product 
BD very rapid. 

Since the matrix G is nonsingular, equation (23) can therefore be written 
as 

^ wheTe C = BD (25) 
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This represents a set of m equations in m unknowns which is exactly equiv- 
alent to the original system of equations. However, if we assume that the known 
vector _f and the solution vector _p can be adequately represented in terms of 
their first n coefficients, then equation (25) can be truncated to give a sys- 
tem of n equations in n unknowns. The precise conditions under which this 
truncation are valid are discussed in detail in section 4.5 at which time an 
improved scheme is developed. However, since this simple truncation is used for 
most of the results presented, further discussion of this point can be postponed 
until section 4.5. 


3.0 APPLICATION OF THEORY TO SMOOTH AIRFOILS 


3.1 Definition of Test Case Airfoils 

For the initial testing of the mode-function method, several special airfoils 
were used in order to give direct control of such factors as: (a) smoothness, 

(b) accuracy of surface definition, (c) exactness of velocity distribution, etc. 
These airfoils were defined by the "pole" theory of James (ref. 5) which uses 
conformal mapping from a unit circle. The two smooth airfoils used are desig- 
nated E33 and El 3 and their shapes are shown in figures 3 and 4. 

The E33 airfoil has a high curvature in the trail ing-edge region but no 
actual discontinuity in surface slope. The source density as a function of the 
arc length will be continuous but there is a sharp peak around the trail ing-edge 
region and so the airfoil does provide a severe test case. The second airfoil, 

£13, does not possess the high curvature of E33 and so it provides a much simpler 
test case. However, comparison of the results for these two airfoils provides a 
useful guide as to how this high curvature affects the rate of convergence of the 
solution for the different methods being tested. It should be noted that, for all 
the results presented here, the point numbering starts from the trailing edge or 
rear stagnation point and proceeds in a counterclockwise direction around the lead- 
ing edge and back to the starting point, (See fig. 2). 

Before presenting the results obtained for these two airfoils it will be 
worthwhile to consider the question of the fittability of the normal velocity 
matrix. As previously noted in section 1.3, the matrix A' in equation (5) is 
not well suited to fitting since it possesses a large spike on the diagonal, but 
the matrix A has a smoother behavior. The 49th row of a 96-element matrix A, 
corresponding to the negative of the normal velocity induced at the 49th control 
point by the other elements (see fig. 3), is shown in figure 5. The Neumann 
matrix used here is based on linear elements with constant source strength on each 
element and so the diagonal term in A is zero giving rise to the spike shown at 
j = 49 in figure 5. This could be removed by including a term dependent on the 
airfoil curvature as in the higher-order method (ref. 6). However, numerical 
results indicate that this effect is of little importance, and it is clearly small 
compared with that found in the A' matrix where the value jumps from about 0.05 
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to 1.0. Therefore, the source distribution results presented here are all based 
on the lower-order matrix formulation. Of greater significance as far as the 
fittability is concerned is the secondary spike caused by the influence of the 
source elements at control points on the opposite surface of the airfoil. This 
can be particularly severe for a thin airfoil where the control point on one 
surface can approach very close to the source elements on the opposite surface. 

Figure 6 illustrates this effect, showing the behavior of the 21st row of the 

matrix A (see fig. 3) for the E33 airfoil. There is no way of removing this 

contribution prior to fitting, and so this effect can be expected to govern the 
number of terms required to obtain an accurate solution. The use of alternative 
singularity mixes might be expected to alleviate this problem since a vorticity 
distribution is known to behave well for thin airfoils. This topic will, however, 
be discussed in greater detail in section 5. 

3.2 Eigenvector Series Results 

It was pointed out in section 2.1 that the number of operations required to 
calculate n eigenvectors associated with a matrix of order m is 0(p n m^) 
where p, the number of iterations required, is about 5 or 6. Therefore, to 
be competitive with a direct solution, the eigenvector method must give acceptable 
answers for n 'v m/(3p) m/15 and to give an improvement in the computer time^ 

a significantly smaller value for n would be required. 

Results obtained for airfoil E33 with 96 points, presented in figures 7-9, 
indicate that this efficiency requirement cannot be achieved. These results 
show the source density and the tangential velocity distribution, plotted against 
the normalized arc length. Thus, for n = 20 or n/m 1/5, the source density 
and the velocity distribution in figures 7 and 8 are inadequate. For n = 30 
the velocity distribution, figure 9, still possesses a significant oscillation. 

The results for the 96-point E13 airfoil are even more unsatisfactory with the source 
density, shown in figures 10 and 11, still oscillating badly even for n = 40, or 
n/m 2/5. 

On the basis of these results, the basic eigenvector method does not offer 
any advantage over the direct matrix solution and, unlike the Fourier expansion 
method, there is no way in which the method can be improved, given the present 
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state of the art. In addition, there does not appear to be any reason 
why the results for three-dimensional problems should be any more favor- 
able. In the light of these results, further work has been concentrated on the 
Fourier method outlined in section 2.2 which does appear to offer significant 
improvements in computation times. The results obtained with the basic Fourier 
method will be shown to be much more encouraging, and in addition, there are 
various techniques which can be applied to improve the convergence and which are 
not available for the eigenvector method. In particular, the use of the FFT 
algorithm enables all of the Fourier coefficients for the matrix to be calculated 

p 

in o(2S m ) operations where S is the sum of all of the prime factors of m, 
m m 

and this knowledge of all the coefficients is crucial to the improved truncation 
scheme to be presented in section 4.5. 

3.3 Sine Series Results 

A good starting point for the sine fit results is the E13 airfoil for which 
the eigenvector series results, figures 10 and 11, were so poor. Thus, figure 12 
shows the 24-term sine series fit to the source density for this airfoil from 
which it is clear that superior results are obtained for a smaller number of terms. 
A similar encouraging trend is observed for the E33 airfoil; figure 13 shows the 
source distribution, plotted against the arc length around the airfoil, obtained 
with a 24- term sine series. Also shown in this figure is the same solution after 
the application of a-smoothing to the expansion coefficients. This technique 
applies an attenuation factor to the Fourier coefficients such that the influence 
of the higher harmonics is decreased. The overall error is increased, but it can 
be seen that the residual oscillation is reduced leading to a more acceptable solu- 
tion. The corresponding tangential velocity distribution for the 24-term expansion 
(without a-smoothing) is shown in figure 14. This truncated solution reproduces 
the sharp peak around the trailing edge remarkably well although there is a sig- 
nificant oscillation present in the distribution. The amplitude of the oscilla- 
tion is considerably reduced by including 32 terms as in figure 15. Alternatively, 
the oscillation could be reduced by the application of a-smoothing, although this 
would reduce the overall agreement. These results, therefore, illustrate the 
potential of the basic sine fitting procedure although the crucial test for the 
technique must be based on realistic airfoil shapes with sharp trailing edges. 
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4.0 APPLICATION TO AIRFOILS WITH SHARP TRAILING EDGES 


4.1 The Basic Sine Series Fit 

As mentioned in the previous section, the combination of the long computation 
times and the poor convergence makes the eigenvector series fit an uneconomic 
proposition. This section will therefore present results obtained by the appli- 
cation of the sine series fit of section 2.2 to airfoils with sharp trailing 
edges. Sections 4.2 to 4.4 will describe some modifications to the basic fitting 
procedure introduced to improve the convergence, and section 4.5 will introduce 
a new matrix truncation scheme which enables an approximate solution to be 
obtained for all of the expansion coefficients of the source density. 

The first sharp airfoil considered is designated BCC#2 and is shown in fig- 
ure 16. This is a very thick airfoil and so it provides a straightforward test 
case with which to continue the validation of the mode-function concept. With 
the freestream parallel to the arrow in figure 16, this airfoil generates no lift 
and at this incidence excellent results were obtained with the basic sine series 
method. The velocity distribution for m = 96 with 8, 12 and 16 terms is shown 
in figures 17-19. A reasonable approximation is obtained with only 8 terms and 
for 16 terms the truncated solution is almost identical to the full solution. 

However, when the airfoil is rotated through about 10®, the results from this 
basic method became unsatisfactory. A sharp oscillation in the velocity distri- 
bution develops near the trailing edge, and this is present even when as many as 
40 terms are used. An examination of the way in which the Kutta condition is 
satisfied for a lifting airfoil reveals the cause of this failure. The singularity 
distribution consists of two components, the first a pure source distribution 
which satisfies the boundary condition on the airfoil due to the onset flow and 
the second a constant vortex plus a varying source distribution which together 
satisfy a zero normal velocity on the airfoil. These components are then combined 
to satisfy the Kutta condition at the trailing edge. For a lifting airfoil, each 
of these components becomes singular at a sharp trailing edge, but in the full 
Neumann solution these two singular terms cancel and the resulting source distribu- 
tion is finite at the trailing edge. However, when the solution is truncated, 
these two singular components are inaccurately approximated and the errors are not 
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eliminated when the combined solution is formed. In order to overcome this 
problem, several modifications to the basic method have been investigated, the 
first two aimed at improving the fitting method in order to cope with singular 
functions while the third method reformulates the Kutta condition so that it is 
satisfied simultaneously with the remaining equations thereby ensuring that the 
source distribution to be fitted is finite at the trailing edge. 

4.2 The Linear Fit 


It will be recalled that the basic sine series fit, equation (18), assumes 
that the function to be fitted is periodic with • This is true in 

theory for a closed body, but if the function to be fitted contains a singular- 
ity at the trailing edge, there will be a large discontinuity between the values 
f^ and f^i+i • However, by subtracting a linear function of the point number, it 
is possible to ensure that the values to be fitted are zero at both ends of the 
interval. Therefore, to fit the vector f_ = (f.| , ... f^^) we first define 

the vector q by 


: , .TLrJ. 

j m - 1 



j = 1 , 2, . . . m 


(26) 


Thus g-j “ Q 3 nd this vector can now be expressed as a sine series, and so 
equation (18) can be replaced by 


Sln[(k-2)(j-l)h] 


(27) 


k=3 


where b^ = f-j , b^ = f^^ and where h is now given by h = 7r/(m — 1). The trans- 
formation matrix g is therefore given by 


G = 


1 

(m-2)/(m-l) 

(m-3)/(m-l) 


0 

l/(m-l) 

2/(m-l) 


;ii 

^21 


l/(m-l) 

0 


(m-2)/(m-l) S 
1 


m-2,1 

0 


where S . . 

* J 


= sin(ijh), 


'1 ,m-2 
’2,m-2 


. S 


m-2,m-2 
0 


( 28 ) 
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This modified fitting procedure is applied both to the influence matrix and the 
column vectors so that the details of the application closely follow those given 

<\j 

in section 2.2 for the basic method. The matrix product D == GG still takes on 
the relatively simple structure 



(29) 


where 

c-| = m(2m-l )/[6(m-l )], 

c^ = m(m-2)/[6(m-l )] 
and 

j cot (j ih), i = 1 , 2, ... m-2 

This can now be substituted into equation (25) to give a nev/ coefficient equation 
which is solved in precisely the same way as in the basic method. 

The results obtained for the lifting airfoil BCC#2 are encouraging with 
the velocity distribution obtained with n = 24 and n = 32 shown in figures 20 
and 21. This represents a significant improvement over the basic method, although 
the oscillation present for the 24-term solution develops rapidly as n is further 
reduced. Comparison of equations (18) and (27) indicates that the number of 
Fourier coefficients has been reduced from (m — 1 ) to (m — 2). However, the 
FFT routine used places some restrictions on this quantity. Therefore, the air- 
foil point number has been increased by one for these and subsequent results. 
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4.3 The Singular Fit 


The breakdown of the results presented above for small values of n indicates 
that the singular behavior of the source density still presents a problem for the 
fitting procedure. It was therefore decided to further modify the technique in 
order to remove the singular contribution before fitting the function. As a very 
crude approximation, the singular behavior of the source density can be expressed 
as the inverse of the point number variable. Reference to equation (26) therefore 
leads to the definition of the vector to be fitted as 




m + 1 - j 



j 1,2, ...m 


(30) 


where b, and b« are defined in terms of f, and f such that g, = g = 0. 
12 1 m ^1 ^m 

Expressing £ in terms of a sine series immediately gives 

m 

^ + I , j b“| -2)(j -l)h] where h = 7:/(m - 1) (31) 

k=3 

Comparison of this with equation (27) shows that in the definition of the trans- 
formation matrix G, equation (28), only the first two columns will be changed. 
Thus 


1/m 

1 

0 . 

0 

1/m-l 

1/2 

$11 . 

• ^1 ,m-2 

1 

1/m 

0 

0 


As with the linear fitting method described in the previous section, this technique 
is applied to each of the vectors and the influence matrix in equation (5) although 
it is primarily aimed at improving the fit to the source density vector. 

It can now be shown that the matrix D = GG will still be given by equation 
(29) with 
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and 




m 




7 




( 33 ) 


The application of this method therefore follows very closely that of the linear 
fitting method once the definition of these constants has been changed. 


This method was applied to the lifting airfoil, BCC#2, defined by 49 points 
and the results are shown in figures 22 and 23 for n = 16 and n = 20. This 
reduction in point number provides a more realistic value for m, but it does 
not lead to a proportionate reduction in n for a given quality solution. The 
res ul ts, therefore, provide a more realistic estimate of the ratio n/m. It can 
be seen that, even for this very crude approximation to the singular behavior of 
the source density, a significant improvement in the truncated solution is 
achieved. A reasonable agreement in the velocity distribution is achieved for 
n = 16 and for n = 20 the results are very close to the full solution except 
near the trailing edge. 


4.4 The Combined Kutta Condition 

The results of the previous two sections illustrate that the basic sine fit- 
ting procedure can be adapted to cope with the singular source distribution which 
occurs in the calculation of the flow over a lifting airfoil. However, figure 23 
still shows a discrepancy between the truncated solution and the full solution 
near the trailing edge. This could be significant if the results were required 
for a boundary-layer calculation and so the precise form of the singularity to 
be removed prior to fitting would need to be investigated further if this method 
were to be pursued. On the other hand, this section considers the possibility of 
reformulating the problem so that the calculated source density is well behaved. 
It was seen in section 4.1 that the basic sine fit gave excellent results for 
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a nonlifting airfoil for which the basic source distribution automatically sat- 
isfies the Kutta condition. Therefore, if the circulation is treated as an 
additional unknown, the Kutta condition can be satisfied simultaneously with the 
source strength equations and the resulting source distribution will be finite 
at the trailing edge. 

This will now give a system of m + 1 equations in m + 1 unknowns which 
could, in principle, be solved by the direct application of any of the fitting 
procedures described above. However, this is not practical since the additional 
values added to the matrix and the column vectors will destroy their continuous 
nature upon which this method relies. For instance, a new column vector would be 
defined consisting of the source strengths and the circulation. Since the cir- 
culation will not be directly related to the adjacent source values, there will 
be a discontinuity in this modified vector. Such discontinuities cannot be ade- 
quately represented by continuous functions and so the total circulation will be 
poorly approximated, and the resulting solution will be slowly convergent in n. 

It is, however, possible to eliminate the circulation numerically and then 
apply the fitting technique to the modified m x m matrix. This approach has 
been applied to the lifting BCC#2 airfoil and the results are shown in figures 
24 and 25 for n = 12 and n - 16. The linear fitting described in section 4.2 
and a-smoothing were used and the resulting velocity distribution is very well 
behaved. Although the agreement near the leading edge is worse than that shown 
in figures 22 and 23, the undesirable oscillation near the trailing edqe has been 
removed. 


4.5 An Improved Matrix Truncation Scheme 

While the results presented so far are encouraging, it should be pointed out 
that these results are all based on the BCC#2 airfoil which, as noted earlier, is 
a very thick airfoil. The next step is therefore to test the method on an airfoil 
with a more realistic thickness, and that used is shown in figure 26. This is 
simply a scaled version of BCC#2, having been derived by dividing the y-ordinate 
of the thick airfoil by 3 and rotating the resulting points through 5®. As the 
thickness of the airfoil is reduced, the peak in the influence matrix resulting 
from the effect of the source elements on the opposite surface, illustrated in 
figure 6, becomes more pronounced. Therefore, this airfoil provides a more severe 
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test case since the fittability of the matrix is impaired. The effect of this can 
be seen in figures 27 and 28 which show the velocity distribution for this airfoil 
using the method of section 4.4 with n = 16 and n = 20, respectively. 

It is clear that even with 20 terms, or n/m 2/5, the results of the trun- 
cated solution are inadequate, while the use of a greater number of terms would 
detract from the computational efficiency of the mode-function method. With the use 
of the FFT algorithm, all of the expansion coefficients for the transformed matrix 
equation (25) are known but the size of the matrix equation to be solved must 
be kept as small as possible in order to reduce the computation time. It was, 
therefore, decided to examine the truncation of the matrix equation with the aim 
of using all of the available information to obtain an improved solution while 
explicitly solving only a small fraction of the equations. 


The transformed coefficient equations to be solved are given by equation 
(25), i.e. 

b = ^ - BD^ (34) 

where D is given either by equation (24) for the basic fit, or by equation 
(29) for the modified fit of sections (4.2) or (4.3). In either case the 
matrix B will consist primarily of Fourier coefficients. Provided that the 
matrix A is sufficiently smooth, the coefficients of the higher harmonics will be 
small compared with those of the lower harmonics. Thus, in general, the most 
significant information contained by the matrix B will be concentrated into the 
first few rows of the first few columns, i.e. the magnitude of the entries of B 
will decrease, in general, as the row or column number increases. Consideration 
of the structure of the matrix D leads to the conclusion that the matrix prod- 
uct C = BD will have the property that, in general, its elements will decrease 
in magnitude as the row number increases. 


In order to express these ideas in mathematical terms, the matrix C and the 
column vectors ^ and ^ in equation (34) can be partitioned to separate the 
first n rows and columns where n < m. Thus 


C = 


'11 



'12 



23 


1 

1 


'^1’ 


and ^ = 


^2 




where C-j-| ts an n x n matrix, while and ^ each have n entries. From 
the preceeding discussion we can assert that, provided n is suitably chosen, 
the largest element in the (m - n) x n matrix C^-j will be small compared with 
the largest element in C-j-j. Similarly the largest element in the (m — n) x (m — n) 
matrix ^22 will also be smal 1 . 

Substitution of equations (35) into equation (34) leads to the two simultane- 
ous matrix equations 


(I C^^)a^ " -1 ^3, 

'"^ 21—1 ^ ^ 22 ^^ ~ —2 

where the symbol I has been used to represent the n x n unit matrix in the 
first equation and the (m — n) x (m — n) unit matrix in the second equation. 
Now if n is "sufficiently large" and ^ and C "sufficiently smooth", then 
^ and C^-j can be neglected and these two equations can be approximated by 


(37) 

^ = 0 

which is the truncated form of the equation which has been used for all the 
results presented so far. At this stage the significance of the dominant 
diagonal becomes clear. If the unit matrices were absent from equation (36) then 
the terms and ^ would all have the same order of magnitude, and the 

approximation ^ = 0 does not necessarily follow. 

Equation (36) is exactly equivalent to the original matrix equation, but 
since 632 is small, it follows that we can make the approximation 

(I - 022)'^ = I + C22 + ... (38) 

and so a^ can be eliminated to give 
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(39) 


[I ^12^^ ^ ^ 1 ” — 1 ^ ^12^^ ^ ^22^^ 

Again, since the matrix C^i is also small, we can introduce the second approx- 
imation 

1 1 
[I ^ ^22^^21^^ "^11^ ^ ^ ^12^^ ^22^^21^^ ”^11^ *‘* 

Therefore, equations (36) can be approximated by 

(I -C^i)a-, = [I + ''22'‘^21^^ -Cn)*’][b^ + Cl2^^ £ 22 )^] 

and 

i2 = (I + C22)(^ + C2^a^) (41) 

The explicit calculation of any of the matrix products involved in the 

right-hand side of equation (41) is out of the question. However, the terms on 

the right-hand side can each be calculated by performing a succession of products 

between matrices and column vectors, together with one matrix inversion required 

to find (I — C-|i)~\ It can be shown therefore that a full solution to equation 

(41) requires 0(n^/3) operations to factorize the matrix (I “C-j-j) together 

with 0(3m^) operations required to perform all of the products necessary to 

solve for a, and The operation count for the solution of the basic trun- 

3 2 

cated equation (37) is 0(n /3), so that for an additional 0(3m ) operations, 
the solution of (41) gives an approximation to all of the expansion coefficients. 

In addition, the solution for the first n coefficients a^^ will also be 
improved since it takes into account more of the available information. 

This scheme has been tested for the thin lifting airfoil, figure 26, and 
excellent results have been obtained. Figures 29 - 31 show the velocity distri- 
bution obtained using n = 8, 12, and 16, respectively, together with the linear 
fit of section 4.2. Thus, even for n = 8, this scheme can produce results which 
could be adequate for many applications, and for n = 16, or n/m 1/3 
the results are very close to the full solution. 

These results serve to demonstrate the potential for this particular formu- 
lation of the mode-function method, although it should be pointed out that they 
do not necessarily represent the optimum. For instance, the only two approximations 
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used are given by equations (38) and (40). However, either of these could be 
made more accurate by the inclusion of higher-order terms, at the expense 
of increasing the number of matrix and vector products which must be calculated. 

A careful balance must be maintained between the need to further reduce the ratio 
of n/m and the additional number of operations incurred in doing so. In addi- 
tion, the results in figures 29- - 31 are based on the linear fitting technique. 
The results of sections 4.3 and 4.4 would indicate that a further improvement 
could be obtained by an improved handling of the singularity which occurs in 
the source distribution. 
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5.0 ALTERNATIVE SINGULARITY MIXES 


5.1 Pure Vortex Formulations 

The results presented in the earlier sections have all been based on the use 
of a source distribution with a constant strength vortex distribution introduced 
in order to provide the circulation about a lifting airfoil. As mentioned in 
section 4.5, the normal-velocity matrix associated with such a singularity distri- 
bution can be difficult to fit for a thin airfoil. However, the matrix does 
possess a dominant diagonal and this has been shown to be important for the suc- 
cess of the improved scheme proposed in section 4.5. On the other hand, the use 
of a pure vortex distribution with an external normal velocity condition is known 
to be well suited to thin airfoil problems. The associated matrix does not 
develop the sharp peaks associated with thin airfoils, but it no longer possesses 
a dominant diagonal. In fact, the values of the elements change sign rapidly 
across the diagonal, and their rate of change increases as the number of ele- 
ments is increased. 

To investigate the relative importance of each of these two factors, the 
basic sine series method together with a pure vortex distribution and external 
normal -velocity condition was applied to the smooth airfoil E33. The vorticity 
distribution used was linearly varying and continuous everywhere except at the 
rear stagnation point, where the first and last values of the vorticity were 
defined to be equal and opposite. 

The velocity distributions derived for m = 96, using 32, 36 and 40 terms, 
are illustrated in figures 32-34. Although the solution, shown in figure 33, 
obtained with 36 terms is good, this agreement must be fortuitous, since the 
additional terms included in figure 34 lead to a deterioration of the solution. 
This is clearly an unsatisfactory situation, and these results should be compared 
with those of figures 14 and 15 which were obtained with the same basic fitting 
method, but using a pure source distribution. The agreement obtained with the 
source method is superior and the behavior of this method was also more predict- 
able, leading to the conclusion that the pure source matrix is more suitable for 
the mode-function approach. In addition, it is not at all clear whether the final 
method described in section 4.5 could be adapted to handle matrices of this form 
and so further work on this particular method was discontinued. 
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An alternative formulation for the pure vortex distribution is provided by 
satisfying the boundary condition that the internal tangential velocity is zero. 

This can be shown to imply that the external normal velocity also vanishes, but 
the numerical formulation is different. In fact, for a vortex density which is 
constant over each element, this problem is numerically identical to the basic 
source formulation used in sections 3 and 4 of this report. However, the applica- 
tion of this boundary condition to the continuous vortex distribution provides a 
convenient way of testing an alternative discretization with the mode-function 
method. The principal difference between the influence matrix for this problem 
and the basic source matrix stems from the choice of the elementary singularities 
used in each case. That selected for the linearly varying vortex distribution 
is taken to be unity at a given point on the airfoil, decreasing linearly to zero 
at the two adjacent points. The singularity is therefore spread out over two 
elements, whereas for the basic source method, the source strength is assumed to 
be unity over the whole length of one single element. This has the effect of 
spreading the “self-influence" spike out over two elements in the matrix thereby 
making it marginally smoother and easier to fit by the mode-function method. On 
the other hand, this self-influence term is no longer given simply by the unit 
matrix and it is no longer feasible to remove it prior to fitting. 

This method has been applied to the E33 airfoil and results for n = 24 and 
n = 32 are shown in figures 35 and 36. Comparison with figures 14 and 15 shows 
roughly the same level of agreement for the basic procedure. As noted above, the 
matrix fitted for this vortex method still contains the unfavorable diagonal 
peak, and so the results indicate that this adverse influence is offset by the 
benefits incurred by spreading the elementary singularities out over two elements. 
However, this particular formulation for the influence matrix renders the applica- 
tion of the improved truncation scheme of section 4.5 much more complicated and 
probably uneconomic. Therefore, this particular formulation has not been pursued 
any further. 

It should be noted, however, that this conclusion does not rule out the use 
of the mode-function method together with a linearly varying singularity distri- 
bution. The problem encountered above is a result of the particular discretization 
selected. Thus, the higher-order method developed by Hess (ref. 6) assumes a linearly 
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varying source distribution, but it allows small discontinuities in the source 
strength from one panel to the next. In this way the elementary singularities 
used are based on single elements and the resulting matrix has essentially the 
same fittability properties as that of the basic method used in section 4.1. 

5.2 Green's Identity Formulations 

The pure source distribution has been shown to be superior to the pure vortex 
distribution as far as the mode -function approach Is concerned. However, for a thin 
airfoil the source strength becomes large and the use of a combined source and 
vortex distribution is attractive. One particularly useful combination is that 
based on Green's identity. This leads to source and vortex densities which are 
equal to the normal and tangential components of the perturbation velocity and the 
result is a source distribution which is more mild than the pure source distribu- 
tion and a vortex distribution which is more mild than the pure vortex distribution. 
This approach has been investigated by Bristow (refs. 7 and 8) who has considered 
two alternative formulations of the problem. The source density follows immedi- 
ately from the known normal velocity and so the problem reduces to that of deter- 
mining the corresponding vortex density. The first approach adopted by Bristow 
(ref. 7) is based on the normal- velocity condition while the second (ref. 8) uses 
an equivalent formulation in which a uniform internal perturbation potential con- 
dition is satisfied. 

The suitability of both of these approaches for use with the mode- function 
procedure has been studied. However, the first method based on the external 
normal velocity condition gives rise to a matrix of the same form as that for the 
pure vortex method considered in section 5.1. This particular formulation was 
found to give unsatisfactory results for the basic fitting procedure and the 
improved technique of section 4.5 cannot be applied, so that it was decided 
not to pursue the corresponding Green's identity approach. 

The second method (ref. 8) with the internal potential condition leads to a 
matrix which is more ammenable to the mode^function approach. The basic sine 
fitting procedure in conjunction with this Green's identity formulation has 
therefore been tested with the BCC#2 airfoil and the results obtained are shown 
in figures 37-40. The first two figures show the results obtained for the 
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nonlifting airfoil with m = 50 for n = 20 and n = 30, respectively, and 
figures 39 and 40 give the results for the lifting airfoil for n = 30 and n = 40. 
The nonlifting results are worse than the corresponding results obtained with the 
pure source distribution for which the basic fitting method was so successful. 

However, it will be recalled that this basic method failed for a lifting air- 
foil (see section 4.1). The Green's identity results do not break down for a lift- 
ing airfoil, but the number of terms required for an accurate solution is clearly 
unacceptable. The problem encountered here is similar to that referred to in 
section 4.4 when the Kutta condition was satisfied simultaneously with the source 
strength equations. In that case- the problem was overcome by numerically eliminat- 
ing the circulation before fitting the matrix. The truncated solution shown in 
figure 39 has the correct overall behavior, apart from close to the trailing edge, 
but the level of the velocity is displaced indicating that the total circulation 
is badly predicted. A similar approach to that adopted in section 4.4 could pre- 
sumably be expected to improve these results also. 

Therefore, this particular form of the Green's identity mix is less suited to 
the mode-function approach than is the pure source distribution. However, the Green's 
identity mix together with a potential boundary condition should, in general, lead 
to a matrix which is diagonally dominant and this should be atmenable to solution 
by the mode -function method. In view of the advantages which this mix can have 
for thin airfoils, further work in this direction would be worthwhile. 
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6.0 APPLICATION TO THREE-DIMENSIONAL PROBLEMS 


6.1 Basic Approach 

The results presented so far are all concerned with two-dimensional applica- 
tions although, as stated in section 1.3, the numerical formulation of a three- 
dimensional Neumann problem is closely related to the two-dimensional case. In 
general, a three-dimensional application will involve a larger number of panels, 
and in addition, the structure of the influence matrix and the singularity vector 
will be rather different. It is, therefore, worthwhile to examine this structure 
before considering a three-dimensional application of the mode-function method. 

A typical three-dimensional configuration is shown in figure 1, and in order 
to formulate the corresponding matrix equation for the singularity strengths, a 
numbering system for the panels has to be selected. A natural choice is to con- 
sider successive streamwise sections so that the wing, for instance, would be 
composed of k strips each consisting of m panels. The corresponding source 
distribution vector would, therefore, consist of k blocks, each containing m 
values associated with the successive streamwise sections. This vector is approx- 
imately periodic in behavior since the values around one strip would be closely 
related to the corresponding values on an adjacent strip. In addition, there will 
be discontinuities in this vector associated with the jump across the trailing 
edge from one strip to the next. In the light of the two-dimensional results, a 
global fit of this function will not produce a very rapidly convergent mode -func- 
tion expansion and a sectional fitting procedure appears to be preferable. 

When applied to the influence matrix, the fitting procedure first fits the 
rows and then the columns of the resulting matrix. Provided that the matrix can 
be held in the core store of the computer, this process can be accomplished effici- 
ently. However, for a three-dimensional problem, this matrix is held on disc 
storage and the penalties involved in accessing this matrix by rows and then by 
columns could be considerable. This also indicates the need to fit the matrix 
in smaller blocks each of which could be held in core and again the use of m x m 
blocks corresponding to the individual streamwise sections appears to be the 
natural choice. 


One final argument against the use of a global fitting method concerns the 
FFT algorithm. The results of section 4.5 clearly demonstrate the need to obtain 
all of the expansion coefficients for both the matrix and the right-hand side, 
and this in turn implies that the FFT algorithm must be used. It has been pointed 
out that this method is able to derive all of the Fourier coefficients for an 
m X m matrix in 0(2 m ) operations where is the sum of all of the prime 
factors of m. Therefore, it is necessary to choose m carefully to avoid includ- 
ing large prime factors which could excessively increase S^. For instance, if 
m itself were a prime number, then S = m and clearly the fitting procedure 
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would be an 0(2m’^) operation and this would be uneconomic compared with the 
direct solutions. In fact, as programmed at Douglas Aircraft Company, the FFT 
routine requires that the value of m (or (m — 1) in the case of the modified fit- 
ting procedures of section 4.2 through 4.5) must be even and a multiple of only 2, 3 
or 5. If the fitting procedure is applied to two-dimensional strips for which 
typically m < 100 this limitation is not unduly restrictive. However, as the 
value of m is increased, the number of acceptable values becomes more sparse. 

For instance, there are only 14 numbers between 2000 and 3000 which would be 
acceptable and this would clearly be far too restrictive. This restriction could 
be partially relaxed by the use of a more general FFT routine but it would still 
be necessary to restrict the size of the maximum prime factor to avoid 
becoming too large, and the use of a global fitting approach may still be too 
restrictive. 

The preferred approach to be adopted for a three-dimensional problem would, 
therefore, be to subdivide the matrix into blocks corresponding to the sectional 
strips. The diagonal blocks will be very closely related to the matrices for the 
two-dimensional problems considered in this report while the off-diagonal blocks 
will define the influence of all the points within one sectional strip upon all 
the control points of another sectional strip. It would not be necessary for these 
individual blocks to have the same dimensions, the only restriction being that 
each of these individual dimensions must be an acceptable value for the FFT 
algorithm. With the unknown source density vector and the given right-hand side 
vector similarly subdivided, the analysis follows very closely that given for 
the two-dimensional method presented here. In particular, the improved scheme 
presented in section 4.5 could be used and so the level of agreement and the 


32 



time saving to be expected from the mode -function approach can be based on this 
two-dimensional study. It will therefore be worthwhile to consider in more detail 
the operation count for a three-dimensional problem handled in this manner. 


6.2 Operation Count for a Three-Dimensional Problem 

Consider the panels to be subdivided into k strips and assume for the sake 

of simplicity that each strip consists of m panels. It has been mentioned that 

the calculation of all of the Fourier coefficients B of an m x m matrix involves 

0(2 m ) operations where is the sum of all of the prime factors of m. 

Typically, for m between 50 and 100, will be between about 11 and 15. Also 

the number of operations required to set up the matrix product C = BO of equa- 

tion (25) will be 0(2m ). Therefore, to find the corresponding matrices for all 

2 ? 

the blocks of the complete matrix will involve 0[2k (S + 1 )m 1 operations. This 
2 2 

will give all km coefficients for the new matrix problem, and so the method of 
section 4.5 can be used. Thus, following section 4.5 let the dimensions of the 
matrix equation to be solved explicitly be kn. Then 0(k^n^/3) operations will 
be required to perform the triangular decomposition of this matrix and a 
further 0(3k m ) operations will be required to complete the solution of equa- 
tions (41). The total number of operations required will therefore be 

2 2 3 3 

0[k (2S + 5)m + k n /3] for the mode-function solution as compared with 

3 3 

0(k m /3) for a direct solution. Thus the ratio between these two quantities 
will be given by 


3 3(2S^ + 5) 


- + 


km 


(42) 


Therefore, this expression gives a guide as to the reduction in computing time 
required to solve the simultaneous equations which could be achieved through the 
use of the mode-function approach. 


The results of section 4.5 indicate that it is realistic to assume that 
n/m = 1/3 although in practice it may be possible to use a smaller value for 
this ratio. For a typical three-dimensional configuration with k = 40 and 
m = 50, = 12 and the value of (42) becomes approximately 1/12. For this 

particular case it would therefore be reasonable to expect the mode-function 
solution to be about 12 times faster than the direct solution. The computing 
effort in this case would be roughly equally divided between the fitting of the 
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matrix and the solution of the resulting equations. As the number of panels is 
further increased, the time required for the fitting becomes a smaller proportion 
of the total time. For such cases it would therefore be worthwhile to improve 
the approximations involved in section 4.5 in order to achieve lower values for 
n/m. Thus for a 4000 panel case with the assumption that n/m * 1/4 the mode- 
function method would be 26 times faster than the direct solution, and an 
8000 panel case would be 38 times faster. On the other hand, an iterative solu- 
tion would require 0(lk m ) operations for one flow condition where I is the 
number of iterations required. For a simple configuration such methods require 
10-15 iterations, although for more complicated configurations such as those 
involving internal flow, more iterations may be required. Thus for one flow 
calculation the iterative method could, at best, be 4-5 times faster than the 
mode- function method for the examples considered above. However, the mode- func- 
tion method does possess several properties of the direct solution which are not 
present in the iterative method. For instance, both the calculation of flows at 
different angles of attack and the simulation of boundary-layer effects require 
at least two flow solutions. Each such solution involves a complete calculation 
for the iterative method whereas for both the direct and the mode- function methods 
the additional solutions can be computed for very little extra cost. In addition, 
the solution time for the iterative method is dependent on the number of itera- 
tions required. As this can vary significantly for different configurations, this 
method is not always predictable. However, provided that n/m is specified before- 
hand, the solution time for the mode=f unction method will be totally predictable. 

Thus, any apparent advantage of the iterative method could very quickly be lost 
in practice. 

The number of operations required to calculate n eigenvectors associated with 

2 

an m X m system of equations is presently 0(pnm ) where p is the number of iter- 
ations required to calculate an eigenvector. Typically p 5 or 6 and so to be 
competitive with the direct solution would require n/m -v 1/15 and to be com- 
petitive with the sine fitting method outlined above would require n/m 1/180. 

In other words, an accurate solution for the 2000 panel example would have to be 
obtained with just 11 eigenvectors. The results of section 3.2 indicate that 
this level of agreement cannot be achieved and so the eigenvector method does not 
provide a viable alternative to the Fourier approach unless a more efficient 
process for calculating the eigenvectors could be developed. 
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7.0 CONCLUSIONS 


The mode-function concept provides a new method for obtaining approximate 
solutions to large systems of linear simultaneous equations. The quantities 
defining many such systems can be considered as discrete values of smoothly 
varying functions and so they can be expressed in terms of a given set of con- 
tinuous mode functions. In this way the original problem can be reformulated 
to give a new set of equations whose solution gives the mode-function expansion 
coefficients for the solution to the original problem. By truncating these 
expansions the number of equations to be solved is reduced and an approximate 
solution to the original set of equations is obtained. 

This report considers the application of this technique to a two-dimensional 
panel method calculation of the inviscid, incompressible flow about an airfoil. 

Two different sets of niode functions are studied, the first are eigenvectors of 
a matrix related to the matrix of the original equations while the second set is 
composed primarily of Fourier functions. An approach to the more time-consuming 
three-dimensional flow problems is proposed. On the basis of the two-dimensional 
results obtained, a substantial reduction relative to the direct solution could be 
achieved in the computer time required to solve the resulting equations. The main 
conclusions to emerge can be summarized as follows: 

1. In their basic form, neither the eigenvector series nor the Fourier series 
approaches to the mode-function concept give adequate results. Thus in 
section 3.2 it was shown that, even with a large number of terms, the 
eigenvector method failed to give adequate agreement for the two smooth air- 
foils considered. In section 4.1, the basic Fourier method was shown to fail 
when applied to a lifting airfoil with a sharp trailing edge. 

2. The computational effort involved in deriving the eigenvectors implies that 
this approach will be an acceptable alternative to the direct solution only 
if a small proportion of the eigenvectors need be calculated. There is no 
way in which the basic theory can be improved to produce a competitive method 
without further analytical work. Therefore, this approach was not pursued 
any further. 
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3. The Fourier series approach is promising since the FFT (Fast Fourier Transform) 
algorithm enables all of the coefficients to be calculated in a time which is 
small compared with a direct matrix solution. This enables modifications to 

be introduced to the basic Fourier method which significantly improve the 
convergence of the solution. Thus the results of section 4.5 show that an 
accurate solution for a realistic airfoil can be obtained by explicitly solv- 
ing less than 1/3 of the original number of equations. 

4. The method of section 4.5 is applicable only to problems whose matrix has a 
dominant diagonal. When applied to the solution of a Neumann problem this 
restriction governs the mathematical and numerical formulation which should 
be used. In particular, an approach based on an unknown vortex distribution 
with an external normal-velocity condition will lead to an integral equation 
of the first kind and this is unsuitable for the mode-function version pre- 
sented here. However, a source distribution formulation leads to an integral 
equation of the second kind and that has been shown to give good results. In 
addition the formulation based on Green's identity together with a potential 
boundary condition leads to a matrix with a similar behavior and the mode- 
function approach will be applicable to such problems. 

5. To avoid fittability problems, the matrix should also be homogeneous in the 
sense that all of the original equations should be of the same form. Thus 
if the Kutta condition or any other equation is to be solved simultaneously 
with the singularity equations then this should be eliminated numerically 
before application of the fitting procedure. 

6. The application of the mode-function method directly to a large matrix 
associated with a three-dimensional problem could involve several numerical 
and computational inefficiencies. These problems could be avoided by 
fitting the matrix in blocks corresponding to two-dimensional sections of 
the configuration. This in turn implies that the results of this two- 
dimensional study should be particularly relevant to the three-dimensional 
problems at which this technique is ultimately aimed. 


7. The operation count presented in Section 6.2 indicates that for a typical 
problem involving 2000 panels the mode-function method could be about 
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12 times faster than the direct solution, and for 4000 panels the method 
could be 26 times faster. Furthermore, this method would possess all of 
the advantages of the direct solution in that the solution time for the 
method would be predictable and additional solutions could be obtained at 
a very small additional cost. 

8. The method, therefore, appears to offer a promising approach to the reduc- 
tion of computing costs for the solution of large three-dimensional Neumann 
problems which are of interest in aircraft design. Alternatively the method 
would offer the capability of handling larger numbers of unknowns for the 
same computational cost. In addition, the method should be applicable to 
the solution of any large system of linear simultaneous equations provided 
only that the matrix has a dominant diagonal, and that it is otherwise 
smooth . 

9. Further work should be undertaken both to consider the application of this 
method to a Green's identity formulation in two-dimensional flow and to 
produce a working three-dimensional method based on either the source 
distribution or the Green's identity method. 
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Figure 1. Representation of a three-dimensional body by surface elements. 
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Figure 2. Linear elements and control points representing a two-dimensional body. 
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Figure 5. The behavior of A. . matrix for i = 49 on the 96-point E33 airfoil 
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